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ABSTRACT 


A  study  of  sedimentation  processes  of  macroscopic,  spherical  particles 
using  a  molecular  dynamics  simulation  was  undertaken.  Calculations  were 
performed  on  particles  of  30,  50  and  100  tin  radius.  For  particle  size 
distributions  where  all  of  the  particles  have  the  same  radius  (monomodal)  it 
was  found  that  the  density  of  the  sediment  is  always  less  than  the  density 
of  an  equivalent  close  packed  structure.  In  addition,  an  increase  in 
particle  radius  leads  to  a  decrease  in  the  sediment  density.  For  bimodal 
particle  size  mixes,  either  30/60,  60/100  or  30/100  radii,  fractionation 
was  found  to  occur.  In  all  cases,  the  larger  particles  settled  before  the 
smaller  particles  giving  rise  to  a  layer  of  small  particles  on  the  top  of  the 
sediment,  and  a  higher  number  ratio  of  larger  particles  in  the  lower 
portions.  . 
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SniULATION  OF  SEDIMENTATION  BY 
MOLECULAR  DYNAMICS 


1.  INTRODUCTION 


Composition  B  is  a  widely  used  military  explosive  that  is  a  60/40  RDX/TNT  mixture. 
T3rpically,  ordnance  is  filled  by  a  casting  process  whereby  Composition  B  is  heated 
above  the  melting  point  of  TNT,  then  poured  into  the  ordnance  shell  and  allowed  to 
cool  and  solidify.  The  Composition  B  mixture  is  typically  heated  10-20 ‘C  above  the 
melting  point  of  TNT  (ca  80 ’C),  i.e.  well  below  the  melting  point  of  RDX  (ca  205*0 
which  decomposes  at  its  melting  point.  At  these  temperatures,  the  solubility  of  RDX 
in  TNT  is  5-6  wt%  (11  hence  this  heated  mixture  is  essentially  a  suspension  of  solid 
RDX  particles  in  liquid  TNT. 

Sedimentation  of  RDX  can  therefore  occur  during  cooling  and 
solidification  of  the  TNT.  Because  batches  of  RDX  have  a  wide  particle  size 
distribution,  fractionation  can  occur  during  sedimentation.  RDX  is  generally  more 
sensitive  to  initiation  than  TNT,  and  therefore  the  fractionation  process  will  produce 
localised  regions  of  higher  RDX  content,  hence  higher  sensitivity.  An  understanding 
of  sedimentation  processes  will  aid  the  production  of  Composition  B  products  and 
other  RDX  based  explosives  where  sedimentation  can  occur,  helping  to  optimise 
homogeneity  of  fills  using  these  explosives. 

In  this  report,  a  computer  simulation  technique  was  used  to  study 
sedimentation  processes.  Since  the  forces  exerted  on  a  particle  in  a  suspension  can 
be  described  readily,  the  method  of  Molecular  Dynamics  (MD)  is  particularly 
suitable.  In  MD  calculations,  the  forces  on  all  particles  are  calculated  and  then  the 
particles  are  moved  according  to  Newton’s  equations  of  motion.  This  process 
develops  from  an  initial  state  of  particles  dispersed  in  a  viscous  medium  subject  to 
gravity  and  the  calculation  is  continued  until  some  final  state  is  reached;  in  this  case 
the  final  state  is  complete  sedimentation.  In  the  calculations,  this  is  assumed 
complete  when  none  of  the  particles  has  significant  motion. 

One  of  the  more  interesting  aspects  of  sedimentation  is  the  effect  of 
particle  size  distribution  on  the  final  density  of  the  sediment.  Information  gained  on 
this  effect  can  be  used  to  tailor  particle  size  distribution  so  as  to  obtain  a  required 
sediment  density. 


2.  MOLECXJLAR  DYNAMICS  (MD)  THEORY 


2.1  Overview 

In  MD  calculations  (see  [2]  for  a  review),  particles  are  moved  according  to  Newton’s 
equations  of  motion  [31,  after  the  forces  on  the  particles  have  been  evaluated.  This 
process  of  force  evaluation/movement  is  repeated  in  cycles  until  the  final  state  is 
achieved.  Consequently,  the  major  features  of  an  MD  code  are; 

(i)  The  force  computation. 

(ii)  Solution  of  Newton’s  equations  of  motion. 

For  computational  convenience  (see  later),  the  force  calculation  and  solution  of  the 
equations  of  motion  are  not  completely  separated.  Some  form  of  "container"  is 
needed  to  describe  the  space  occupied  by  the  particles.  For  sedimentation,  a  square 
beaker  has  been  used  where  the  X  and  Y  directions  correspond  to  the  width  and  depth 
of  the  container.  The  Z  direction  is  aligned  with  the  direction  of  gravity  and 
corresponds  to  the  height  of  the  container.  In  order  to  make  this  beaker  behave  like  a 
real  beaker,  certain  constraints  are  placed  on  the  X,Y  and  Z  directions  and  these 
constraints  will  be  described  in  detail  in  section  2.2. 


2.2  Faroe  Calculation 

'The  important  forces  in  a  sedimentation  process  are  viscosity,  gravity,  interparticle 
forces,  particle /floor  forces  and  buoyancy. 


The  magnitude  of  the  viscous  force  on  a  spherical  particle  is  given  by  l41; 


F  =  -Girnr.v.  (1) 

where  n  is  the  liquid  phase  viscosity,  rj  is  the  radius  of  particle  i  and  V;  is  the  velocity 
of  particle  i.  The  viscous  force  acts  in  the  direction  opposite  to  the  direction  of  the 
motion  of  the  particle. 


Gravitational  force  is  given  by: 


m.g 


(2) 


where  mj  is  the  mass  of  the  i’th  particle  and  g  is  the  gravitational  field  strength. 

In  the  present  calculations,  the  gravitational  force  is  defined  as  being  in  the  -Z 
direction  (i.e.  g=  I0,0,g^]). 

'The  interparticle  forces  describe  the  interaction  between  particles.  Here, 
the  particles  are  assumed  to  be  spherical  "soft"  spheres.  A  "soft"  sphere  model  is 
necessary  to  overcome  numerical  stability  problems  as  well  as  reducing  particle 
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overlap.  There  is  no  long-range  attraction  (e.g.  opposite  charge  attraction)  or 
cohesion /adhesion  force.  A  truncated  Lennard-Jones  12-6  type  potential  proved  to  be 
the  most  suitable  and  is  given  by: 


F. 


ij 


=  4...( 


_ LL 

13 
r .  . 


606. 
_ u 


)  d 


(3) 


where  n :  is  the  interparticle  distance,  e .  ,  and  a.  .  are  constants  for  pair  (i,j)  and  d  is 
a  unit  vector  along  the  interparticle  line*.'^  The  direction  of  this  force  is  along  the 
interparticle  vector.  The  range  of  Fj  j  is  0  to  r^  +  rj,  that  is  from  zero  to  the 
sum  of  the  particle  radii,  and  for  r.  .  >  r.  +  r  .  we  assume  Fjj  =  0.  This  means  the 
force  starts  to  operate  when  the  parlicles^starf'to  touch.  SeeTigure  1  for  a  typical 
plot  of  this  function. 


In  the  simulation,  there  is  a  restriction  on  the  number  of  particles  that  can 
be  included  because  of  the  limitations  imposed  by  the  storage  capabilities  of  the 
computer  being  used.  Because  of  this  limit,  a  technique  called  minimum  imag^g  is 
employed  to  overcome  undesirable  edge  effects.  With  minimum  imaging,  the  box  of 
particles  being  studied  is  surrounded  by  mirror  images  of  itself  in  the  X  and  Y 
directions.  When  interparticle  forces  are  calculated  between  a  pair  (i,j),  the  i’th 
particle  is  deemed  to  be  in  the  original  box  and  the  j’th  particle  is  then  chosen  from 
the  original  or  images  such  that  the  distance  is  a  minimum.  Another  way  of 
describing  this  is  that  a  particle  leaving  the  box  in,  say,  the  +X  direction  is  treated  as 
entering  the  box  from  the  -X  direction.  Thus  minimum  imaging  also  describes 
"periodic  boundary  conditions". 


A  particle-floor  potential  is  needed  between  the  particles  and  the  bottom 
of  the  container.  The  form  of  this  potential  was  taken  as: 


j.floor^^ 


(4) 


where  r .  is  the  distance  between  particle  i  and  the  bottom  of  the  container,  K;  is 
a  constant  for  particle  type  i  and  k  is  a  unit  vector  in  the  Z  direction.  The  form  of 
this  potential  is  basically  the  repulsive  part  of  equation  (3)  and  was  chosen  for 
consistency  as  well  as  its  mjccess  in  providing  an  adequate  description  of  the  bottom 
boun^n^.^  The  range  of  F/  _  j  i^jp  to  rj,  i.e.  from  zero  to  the  particle  radius,  and 
forr.  >  we  assume  F:*  =0.  Figure  2  gives  a  typical  representation  of 

this  function,  it  should  be  noied  that  minimum  imaging  is  not  employed  in  the 
direction  of  gravity  (i.e.  the  Z  direction). 


Strictly  speaking,  a  buoyancy  force  should  also  be  incorporated  into  the 
calculations.  However,  at  this  stage  we  are  primarily  interested  in  the  way  particles 
pack  during  sedimentation  and  less  interested  in  rates  of  sedimentation.  Including 
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buoyancy  will  only  reduce  the  effect  of  gravity,  since  it  operates  in  the  exact  opposite 
direction,  and  is  not  expected  to  affect  the  mode  of  packing.  Consequently,  buoyancy 
has  been  omitted  from  the  calculations  presented  here. 


The  effect  of  Brownian  motion  on  the  sedimenting  particles  has  not  been 
included  in  the  calculation  as  the  particle  sizes  considered  here  are  large  enough  for  it 
to  be  neglected.  A  simple  criterion  for  the  neglect  of  Brownian  motion  is  that  the 
particle  Peclet  number  be  much  larger  than  unity  [51.  For  the  system  considered  here 
this  will  be  the  case  provided  the  particle  number  is  greater  than  1  fjm.  The  smallest 
particle  radius  considered  here  is  30  fxn,  hence  Brownian  motion  has  not  been  included 
in  the  calculation. 


The  total  force  on  a  particle  is  therefore  given  by: 


F.  =  m.a.  =  I  F.  .  +  m.g  -  Bnqr.v.  +  (5) 

1  1  j  j=i  *  'ill 

j-i 

where  N  is  the  number  of  particles,  Vj  is  the  velocity  and  blj  is  the  acceleration  of 
particle  i. 


2.3  Solving  the  equations  of  motion 

For  a  molecular  dynamics  simulation,  the  velocity  and  position  at  one  time  step  are 
related  to  the  velocities  and  positions  at  other  time  steps  via  the  Verlet  equations  161: 

+  At)  ~  r.(t  -  At))  (6) 

1  2At  1  1 

Tjit  +  At)  =  2r.(t)  -  r.(t  -  At)  +  aj(t)At^  (7) 


where  vj(t)  is  the  velocity  vector  of  particle  i  at  time  t,  aj(t)  is  the  acceleration 
vector  of  particle  i  at  time  t  and  rj(t)  is  the  position  vector  of  particle  i  at  time  t.  To 
calculate  r.  (t  +  At),  a.(t)  must  be  known.  From  equation  (6),  for  aj(t)  we  need 
v^(t)  which,*in  turn,  is  de^ndent  on  rj(t  +  At) .  This  difficulty  can  be  overcome  as 
follows  (considering  only  the  x-component  for  one  particle).  Equation  (5)  can  be 
written  as: 


T  F.  .  +  m.g  +  Ft 

iix  i“x 


floor 


Srinr  .  dx 
1  _ n 

m.  dt 
1 


(8) 


where  Xjj  is  the  x-coordinate  of  particle  i  at  time  t.  Since  Fj;,  and  F^  are  only 
dependent  on  the  particle  positions  at  time  t  (which  are  known)  and  m^*^  is  a  constant, 
equation  (8)  then  becomes: 
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A  - 


dx 

n 

“  dt 


(9) 


where 


A 


—  (  I  F.  .  +  m.g  +  f!^“°’‘) 

m.  ^  ijx  i**x  IX 
i 


a 


Grrrjr. 

m. 

1 


(10) 


(11) 


E)quation  (7)  can  be  written  as  (x-component  only): 


d  X 


n+1 


=  2x  -  X  ,  + 
n  n-1 


dt^ 


(12) 


where  is  the  particle  x-coordinate  at  time  t  +  At  and  x^-l particle 
x-coordinate  at  time  =  t  -  At .  Equation  (6)  can  then  be  simplified  using  equations  (9) 
and  (12): 


dt 


■TTT  (x  1  -  X  -1  ) 
2At  n-t-l  n-1 


1  *n  2 

■577  (2x  -  2x  ,  +  - n  At  ) 

2At  n  n-1  ^^2 

^n  -  ^-1  At  ila 
At  2  ^^2 


^n  ~  ^n-1  M  A  _  “M.  ^ 
At  ^2  2  dt 


(13) 


Collecting  dxjj/dt  on  the  RHS  gives: 


dt 


_ 1 _ 

(1  +  a(At/2)) 


-  Vl 

At 


(14) 


Equation  (14)  can  be  evaluated  given  Xjj  and  Xj._j^  only.  Initially,  particle  positions  and 
velocities  are  assigned  then  Xj^  is  approximated  via  : 
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(15) 


Subsequent  MD  cycles  follow  the  sequence  below: 

1/  Evaluate  A  using  equation  (10) 

2/  Solve  equation  (14)  giving  the  particle  velocities  dx  ^t 
3/  Solve  equation  (9)  giving  the  particle  accelerations  a^x^/dt^ 
4/  Solve  equation  (12)  giving  Xj,^j 
5/  Increase  the  time  step  by  an  increment 
6/  Start  again  at  step  1. 

In  this  scheme,  the  force  evaluation  and  solution  of  the  equations  of  motion 
are  interleaved  to  overcome  the  velocity  dependence  of  the  force. 


One  further  point  to  note  is  the  use  of  periodic  boundary  conditions  in  the 
X  and  Y  directions,  as  discussed  in  Section  2.2.  If  a  particle  moves  outside 
of  the  X  or  Y  limits  of  the  container  it  is  removed  and  replaced  by  an  identical 
particle  (with  the  same  velocity  vector)  entering  the  container  from  the  opposite  side. 


3.  VOLUME  CALCULATIONS 


In  order  to  compare  results  for  different  particle  size  distributions, 
a  method  for  calculating  the  sediment  volume  is  required.  Given  the  box  dimensions, 
the  particle  coordinates  and  radii,  the  volume  can  be  calculated  by  evaluating  the 
volume  under  the  top  surface  of  the  particles.  The  top  surface  can  be  defined  by 
lowering  a  grid  of  small  squares  onto  the  sediment  until  all  of  the  squares  have 
intersected  a  particle  (Figure  3).  The  volume  is  then  the  sum  of  the  volumes  under 
every  square,  and  shall  be  called  Vwq  for  volume  from  molecular  dynamics 
calculation.  For  comparison,  we  define  two  other  volumes: 


V 

8 


477  ^  3 

3“  ^  ''i 

i=l  ' 


V 

c 


N 

4V2  I 
i=l 


3 


r . 


1 


(16) 


(17) 


Vg  corresponds  to  the  total  volume  of  the  spheres  only  and  ignores  any  interstices 
created  when  spheres  pack  together.  is  the  close  packed  volume  for  either  HOP 
(hexagonal  close  packing)  or  FCC  (face  centered  cubic).  It  is  expected  that: 
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2  V 


> 


V 
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so  Vp  and  Vg  provide  a  check  on  the  value.  If  is  equal  to  Vp  then  the 
particles  must  be  close  packed.  The  use  of  Vp  is  only  meaningful  for  situations  where 
all  the  particles  are  the  same  size. 


The  density  of  the  sediment,  pj^,  is  given  by: 


m. 

1=1  1 


WD 


(4rrp  )  t:  ,r. 
O  1=1  1 


HD 


V 


s 


(18) 


where  m^  is  the  mass  of  the  i’th  particle  and  p  is  the  particle  material  (crystal) 
density.  ° 


When  comparing  results  wnere  the  only  difference  is  the  particle  size 
distribution,  we  can  define  a  relative  density: 


P 


r 


(19) 


must  be  s  1,  and  in  general  will  be  significantly  less  than  unity.  The  nearer  to 
unity  it  becomes,  the  more  closely  packed  is  the  sediment. 


4.  RESULTS 


4.1  Computational  considerations 

One  of  the  less  desirable  features  of  a  molecular  dynamics  code  is  its  order  N  usage  of 
CPU  time  where  N  is  the  number  of  particles  in  the  simulation.  This  means  that  an 
increase  in  the  number  of  particles  by  a  factor  of  two  will  give  a  factor  of  four 
increase  in  the  CPU  ♦ime  needed  to  perform  the  simulation.  It  is  possible  to  write 
MD  codes  that  give  improved  performance  (71  but  these  are  more  complicated  and 
take  longer  to  develop.  Since  it  was  not  clear  if  order  N  performance  would  be 
prohibitive,  the  extra  effort  to  write  a  more  efficient  code  was  not  undertaken. 
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Apart  from  changing  the  basic  algorithm  used  in  the  MD  code, 
improvements  in  time  efficiency  can  be  obtained  by  implementing  other  innovations. 
In  this  work,  the  interparticle  and  particle/floor  forces  are  tabulated.  This  removes 
the  time  expensive  process  of  r^^  and/or  r^^  evaluation  every  time  the  force  is 
calculated  and  replaces  it  with  a  look-up  table.  Interpolation  is  not  necessary  due  to 
the  high  resolution  of  the  look-up  table.  The  At  used  in  Equations  (6)  and  (7)  can  be 
increased  to  reduce  the  number  of  cycles  needed  for  complete  sedimentation. 
However,  if  At  is  increased  too  much  the  calculation  can  become  unstable.  This 
instability  arises  from  excessive  particle  motion  in  a  given  cycle.  The  large  motions 
lead  to  excessive  particle  overlap  which  in  turn  cause  huge  forces  on  the  overlapping 
particles  that  finally  leads  to  particle  movements  many  times  the  size  of  the  box 
dimensions.  Such  movements  cannot  be  handled  adequately  and  lead  to  numerical 
problems.  Consequently,  At  must  be  chosen  very  carefully  to  avoid  instabilities  yet 
allow  as  lOw  cycles  as  possible  for  complete  sedimentation.  The  severity  of  this 
instability  problem  can  also  be  influenced  by  the  choice  of  interparticle  potential. 


Initially,  the  particles  are  placed  randomly  within  a  box  and  given  a 
random  velocity.  To  obtain  statistically  meaningful  results,  calculations  should  be 
performed  a  number  of  times  with  different  initial  conditions.  The  pseudo  random 
number  generator  used  to  provide  these  initial  conditions  is  based  on  the  recurrence 
relation  [81; 


p(i)  =  AfKi  -  Dmod  M 


(20) 


where  p(i)  is  the  i’th  pseudo  random  number  generated,  M  is  the  period  of  the  pseudo 
random  number  generator  and  A  is  a  constant.  The  p(i)  are  distributed  between  0  and 
M.  Therefore  to  obtain  random  numbers  between  0  and  1  all  that  is  required  is 
division  by  M. 

In  addition,  instead  of  using  the  viscosity  of  liquid  TNT  at  an  appropriate 
temperature,  the  viscosity  of  water  at  25' C  (17=  8.904  x  10'^  Pa.s  191)  was  used 
throughout  the  calculations.  The  reason  for  this  comes  back  to  computer  time;  TNT 
has  a  higher  viscosity  so  more  cycles  are  needed  for  complete  sedimentation  and 
hence  more  CPU  time  wilLbe  used.  The  particle  material  density  was  also  given  an 
arbitrary  value  of  2  Mg/m".  This  value  doesn’t  reflect  the  density  of  any  material  in 
particular,  but  is  a  value  that  has  the  same  order  of  magnitude  as  the  density  of 
typical  explosive  solids. 


Values  of  e  and  a  can  be  obtained  in  the  literature  [101  for  atoms  and 
molecules.  However,  e  and  o  values  for  macroscopic  particles  like  those  used  in  this 
study  are  not  readily  available.  The  values  used  here  have  been  chosen  somewhat 
arbitrarily  but  satisfy  the  required  conditions  fairly  well,  i.e.  the  particles  behave  as 
"soft"  spheres.  K  values  are  also  assigned  in  a  similar  fashion.  Tables  1,  2  and  3  list 
values  of  e ,  o  and  K  used  in  this  work.  The  results  we  obtained  were  found  not  to  be 
sensitive  to  the  choice  of  e ,  a  or  K. 


4.2  All  particles  the  same  size  (monomodal  blend) 

Initial  calculations  were  performed  for  suspensions  where  all  of  the  particles  have  the 
same  radius.  These  results  provide  an  indicator  as  to  the  reliability  of  the  code  as 
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well  as  providing  a  simple  starting  point  for  more  involved  calculations.  Three 
particle  sizes  were  investigated;  r  =  30,  50  and  100  pm  (Tables  4,  5  and  6).  In  order 
to  obtain  statistically  meaningful  results,  similar  calculations  were  performed  a 
number  of  times  for  each  radius.  The  only  difference  between  these  calculations  is 
the  starting  positions  and  velocities  of  the  particles.  Insp>ection  of  the  p  values 
(Table  7)  for  the  three  radii  show  that  p  tends  to  decrease  slightly  with  increasing 
radius.  This  suggests  that  smaller  partKilcs  tend  to  pack  more  closely  than  larger 
particles. 


The  size  of  the  calculation,  i.e.  the  number  of  particles,  can  have  a  marked 
affect  on  the  results  of  MD  calculations  especially  if  periodic  boundaries  are  not 
used.  The  size  effect  was  investigated  by  repeating  the  50  }Jxi  radius  calculations 
with  N  =  216  (Table  8).  There  is  no  statistically  significant  difference  in  the  Vwjj 
value  for  N  =  125  and  N  =  216.  This  indicates  that  N  =  125  is  sufficient  to  obtam 
reliable  results  for  monomodal  particle  size  distributioiL  Comparison  of  the  average 
times  for  N  =  125  and  N  =  216  demonstrate  the  order  performance,  (216/125)^  ...  3; 
one  would  therefore  expect  a  factor  of  three  increase  in  CPU  time  for  the  N  =  216  run 
as  compared  to  the  N  =  125  run.  From  Table  5,  the  N  =  125  run  takes  about 
42  minutes  and  from  Table  8  the  N  =  216  run  takes  about  125  minutes  so  the  expected 
time  increase  is  the  same  as  the  observed  time  increase. 


4.3  Particles  of  two  sizes  (bimodal  blends) 

Three  mixes  were  investigated  here,  namely,  30/60,  60/100  and  30/100  ^n  mixes.  The 
number  ratio  of  the  small  to  large  particles  was  held  fixed  at  6:3.  The  results  are 
presented  in  Tables  9,  10  and  11.  With  a  mix,  there  is  the  possibility  of  fractionation 
occurring  during  sedimentation  therefore  plots  of  the  probability  of  locating  a  particle 
centre  as  a  function  of  height  (in  the  Z  direction)  were  obtained.*  To  improve  the 
statistics  of  these  plots,  the  final  particle  coordinates  of  all  of  the  runs  of  a  particular 
mix  were  combined  and  then  analysed.  The  plots  obtained  are  shown  in  Figures  4,  5 
and  6. 


The  bimodal  calculations  were  more  susceptible  to  instability  problems 
than  the  monomodal  of  section  4.2,  requiring  a  smaller  it  in  conjunction  with  more 
cycles  for  the  simulations.  In  addition,  the  greater  the  difference  in  the  particle  size, 
the  greater  was  this  instability.  One  possible  explanation  for  this  lies  in  the  relative 
terminal  velocities  of  the  particles.  By  considering  the  free  fall  of  a  particle  in  the 
suspension  (free  fall  here  means  there  are  no  interparticle  or  particle /floor  forces  on 
the  particle)  it  can  be  shown  that  the  terminal  velocity  of  a  particle  is  proportional  to 
the  square  of  the  particle  radius.  Thus,  larger  particles  move  with  a  higher  velocity 
and  will  move  greater  distances  in  a  given  cycle  than  smaller  particles.  It  would  then 
seem  possible,  within  a  single  time  step,  for  a  large  particle  to  considerably  overlap  a 
smaller  particle. 


*  In  fact,  sedimentation  is  a  well-known  technique  for  particle  sizing 
measurements  [111. 
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6.  DISCUSSION 


The  most  important  parameter  calculated  for  this  study  was  Vwjj,  and  it  was 
important  to  ensure  that  enough  independent  estimates  were  obtained  to  provide  a 
representative  sample.  To  gauge  this,  a  t-test  was  performed  for  the  values  that 
gave  the  largest  standard  deviation,  that  is  the  100  (Jn  runs  (see  Table  6).  The  six 
values  were  split  into  two  groups  of  three,  one  group  containing  the  three  highest 
values  and  the  other  group  with  the  three  lowest.  The  test  hypothesis  that  the 
difference  in  the  average  value  of  the  two  groups  is  zero  results  in  t  becoming  equal 
to  2.55,  thus  the  hypothesis  is  accepted  at  the  99.5%  confidence  limit.  'Diis  worst 
case  test  suggests  that  six  independent  estimates  provide  a  sufficient  coverage  of  the 
range  of  possible  values. 


In  section  4.2  (see  also  Table  7)  it  was  noted  that  smaller  particles  tend  to 
pack  more  closely  than  larger  particles.  If  the  spheres  were  close  packed  then  p  = 
0.74  irrespective  of  particle  size  so  it  seems  reasonable  to  assume  that  p  should^e 
independent  of  particle  size  (if  all  of  the  particles  have  the  same  radius  in  a  given  run) 
for  a  given  mode  of  packing.  The  results  in  Table  7  suggest  that  the  30  tm  particles 
pack  in  a  different  way  to  the  100  tm  particles  and  this  difference  can  be  explained  by 
the  r^  dependence  of  the  terminal  velocity.  The  smaller  particles  have  a  slower 
terminal  velocity  so  they  have  more  time  to  rearrange  and  adopt  a  more  efficient 
packing  structure  than  the  larger  particles.  Alternatively,  the  difference  could  be 
due  to  the  interparticle  potential  that  might  cause  the  interaction  between  particles 
to  be  radius  dependent  and  thus  affect  the  mode  of  packing.  However,  this  is  less 
likely  since  the  form  of  this  function  and  the  method  for  choosing  the  parameters  is 
the  same,  irrespective  of  particle  radius. 


The  possibility  of  fractionation  in  the  bimodal  simulations  is  evident  from 
the  tr  dependence  of  the  terminal  velocity.  The  larger  particles  tend  to  settle  before 
the  smaller  particles  because  the  larger  particles  have  a  higher  terminal  velocity. 
Plots  of  probability  as  a  function  of  height  for  the  mix  simulations  show  the  extent  of 
fractionation.  Elach  plot  has  a  separate  curve  for  every  distinct  particle  size  in  the 
mix,  therefore  for  the  30/60  ^^n  simulation  (Figure  4)  there  is  one  curve  for  the 
30  Mn  particles  and  one  for  the  60  ^n  particles.  The  same  applies  for  the 
50/100  fxn  (Figure  6)  and  the  30/100  fxn  (Figure  6)  simulations. 

A  common  feature  of  these  plots  is  a  pair  of  well  defined,  sharp  peaks  at 
low  heights.  For  example,  the  30/60  plot  (Figure  4)  has  a  sharp  peak  at  26  ^n  for 
the  30  tin  particles  and  another  at  46  ^n  for  the  60  ;xn  particles,  lliese  peaks 
correspond  to  particles  in  contact  with  the  floor  of  the  container.  The  fact  that 
these  peaks  occur  at  heights  slightly  lower  than  the  corresponding  particle  radius 
indicate  a  small  amount  of  overlap  between  the  floor  and  the  particles  (i.e.  squashing). 

The  tendency  for  the  smaller  particles  to  settle  after  the  larger  particles 
is  shown  by  the  plots.  For  the  30/60  ^fn  simulation  (Figure  4)  the  30  tin  probability 
curve  extends  out  to  360  jxn  whereas  the  60  ^n  curve  only  goes  out  to  310  ^In.  Thus, 
beyond  310  «4n  there  is  zero  probability  of  finding  a  50  particle  and  only 
30  fin  particles  will  be  found  in  the  range  310  tm  to  360  fm,  i.e.  a  layer  of 
30  im  particles  will  be  formed  on  top  of  the  sediment.  Similarly,  for  the 
60/100  fin  simulation  (Figure  6)  100  fin  particles  are  not  found  above  678  fin  and  only 
50  fm  particles  are  found  between  678  fin  and  623  fin.  The  30/100  fin  simulation 
(Figure  6)  is  not  as  clear  cut  since  the  probability  of  finding  a  100  fin  particle  above 
470  fm  is  not  zero  but  since  it  is  very  small  we  can  take  470  as  the  upper  limit  for  the 
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100  fin  particles  here.  Only  30  fin  particles  will  be  found  between  470  fin  and 
530  fin. 


6.  CONCLUSION 


For  a  particle  size  distribution  where  the  particles  all  have  the  same  radius,  the 
simulations  indicate  that  the  packing  density  is  always  less  than  the  equivalent  close 
packed  structure.  For  a  close  packed  structure  p  =  0.74  whereas  the  simulations 
gave  p  =  0.72,  0.70,  0.68  for  particle  radius  of  36,  50,  100  fin  respectively.  Thus, 
the  packing  density  was  found  to  be  particle  size  dependent;  an  increase  in  particle 
radius  gave  a  decrease  in  sediment  density.  When  the  mix  consists  of  particles  of  two 
distinct  sizes,  fractionation  was  found  to  occur.  The  larger  particles  have  a  higher 
terminal  velocity  than  the  smaller  particles  so  they  tend  to  settle  before  the  smaller 
particles.  This  fractionation  gives  rise  to  a  layer  of  smaller  particles  on  the  top  of 
the  sediment.  This  type  of  inhomogeneity  has  important  consequences  for  the 
sensitivity  and  mechanical  properties  of  an  explosive  where  sedimentation  may  have 
occurred  during  production. 

Understanding  packing  density  for  explosive  fillers  is  important  for 
explosive  formulation,  production  and  filling.  TTie  scope  of  the  calculations  presented 
in  this  paper  has  been  limited  by  the  computer  time  needed  to  perform  the 
calculations.  For  example,  approximately  22  hours  of  VAX  8700  CPU  time  is  required 
to  complete  six  runs  of  200  particles  (mixed  sizes).  In  order  to  perform  these 
calculations  routinely  and  study  larger  systems  (ie  more  complex  particle  size 
distributions  and  more  particles),  improvements  in  the  time  efficiency  of  the  code  will 
be  essential.  It  is  hoped  that  incorporating  a  vectorizable  near-neighbours  algorithm 
[61  into  the  code  will  allow  a  comprehensive  series  of  calculations  to  be  performed 
without  too  great  a  demand  on  computer  resources. 
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A 


p(i) 

M 


total  volume  of  particles 
close  packed  volume  of  spheres 
volume  of  sediment  from  MD  simulation 
density  of  sediment  from  MD  simulation 
relative  density  of  sediment  to  particle  material 
A  is  a  constant 

p(i)  is  the  i’th  pseudo  random  number  generated 
M  is  the  period  of  the  pseudo  random  number  generator 
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Table  1 


e  values  used  in  the  interparticle  force  calculation  for  particles  of  radius 
30,  50  and  100  tia. 


Table  2 


e 

X  10^^  (J) 

30 

50 

100 

30 

0.473 

1.50 

10.4 

50 

1.50 

3.65 

18.5 

100 

10.4 

18.5 

58.4 

o  values  used  in  the  interparticle  force  calculation  for  particles  of  radius 
30,  50  and  100  ^In. 


a  X  10®  (m) 


(m) 

30 

50 

100 

30 

14.7 

19.6 

31.9 

50 

19.6 

24.5 

36.8 

100 

31.9 

36.8 

49.0 

Table  3  K  values  used  in  the  particle /floor  force  calculation 


Radius 

(jm) 

30 

60 

100 

K 

3.64  X  lO"®® 

1.25  X  lO"®"* 

8.22  X  lO"®® 

21 


Table  4  Results  for  30  fin  particles  simulation 
Time  step  =  250  fiB 
Cycles  =  6000 

XY  dimensions  =  300  tm  by  300  fin 
Starting  Z  =  417  tm 
Particle  radius  =  30  jm 
No.  of  particles  =  125 


Run  No. 

Vmd 

(X  10'^^  m^) 

Run  time 

(mins) 

1 

1.969 

46 

2 

1.976 

43 

3 

1.970 

42 

4 

1.970 

42 

5 

1.978 

42 

6 

1.970 

46 

Average 

1.972 

43 

Std  Dev 

0.004 

2 

Table  5  Results  for  50  fm  particles  simulation 
Time  step  °  250  fis 
Cycles  =  6000 

XY  dimensions  e  600  fin  by  500  im 
Starting  Z  =  694  fm 
Particle  radius  =  50  fjn 
No.  of  particles  =  125 


Run  No. 

'^MD 

(X  10"^^  m®) 

Run  time 

(mins) 

1 

9.41 

45 

2 

9.49 

43 

3 

9.38 

41 

4 

9.00 

41 

5 

9.21 

42 

6 

9.60 

42 

7 

9.08 

40 

Average 

9.31 

42 

Std  Dev 

0.2 

2 

22 


j 


Table  7 


Table  6  Results  for  100  ixa  particles  simulation 
Time  step  =  260  jjb 
Cycles  =  6000 

XY  dimensions  =  1000  fin  by  1000 
Starting  Z  =  1389  isn 
Particle  radius  =  100  pea 
No.  of  particles  =  125 


Run  No. 

(X  10"^^  m®) 

Run  time 

(mins) 

1 

73.71 

41 

2 

80.51 

41 

3 

79.77 

41 

4 

79.31 

41 

5 

78.66 

40 

6 

75.98 

40 

Average 

78 

41 

Std  Dev 

2 

0.5 

values  for  the  one-particle  size  distribution  simulations 


Radius 

(pen) 

N 

Pr 

30 

126 

0.72 

60 

126 

0.70 

60 

216 

0.70 

100 

126 

0.67 

23 


Table  8 


Results  for  50  ^fn  particles  simulation  with  an  increased 
number  of  particles 
Time  step  =  250  ps 
Cycles  =  6000 

XY  dimensions  =  600  ^n  by  600  iJn 
Starting  Z  =  882  tjoa 
Particle  radius  =  50  tia 
No.  of  particles  =  216 


Run  No. 

'^MD 

(X  10"^^  m^) 

Run  time 

(mins) 

1 

15.80 

117 

2 

16.80 

127 

3 

16.41 

128 

4 

16.43 

128 

5 

16.33 

118 

6 

16.41 

129 

Average 

16.2 

125 

Std  Dev 

0.3 

5 

Table  9  Results  for  30/50  tm  particles  mix  simulation 
Time  step  =  200  ps 
Cycles  =  12000 

XY  dimensions  =  500  nn  by  500  im 
Starting  Z  =  850  tm 

Particle  radius  =  30  jm  No.  of  particles  =  125 
Particle  radius  =  50  Mn  No.  of  particles  =  75 


Run  No. 

Vmd 

(X  lO'^l  m^) 

Run  time 

(mins) 

1 

7.639 

240 

2 

7.661 

210 

3 

7.691 

202 

4 

7.667 

206 

5 

7.663 

203 

6 

7.678 

201 

Average 

7.66 

210 

Std  Dev 


0.03 
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Table  10:  Results  for  30/100  tin  particles  mix  simulation 
Time  step  ==  150  fjs 
Cycles  =  14000 

XY  dimensions  =  1000  im  by  1000  im 
Starting  Z  =  1306  im 

Particle  radius  =  30  tm  No.  of  particles  =  125 
Particle  radius  =  100  tm  No.  of  particles  =  75 


Run  No. 

'^MD 

(x  10"^^  m^) 

Run  time 

(mins) 

1 

45.85 

240 

2 

45.43 

210 

3 

48.05 

202 

4 

46.40 

206 

5 

48.44 

203 

6 

48.48 

201 

Average 

47 

237 

Std  Dev 

1 

9 

Table  11:  Results  for  50/100  tm  particles  mix  simulation 
Time  step  =  200  ps 
Cycles  =  12000 

XY  dimensions  =  1000  tm  by  1000  tin 
Starting  Z  =  1510  tm 

Particle  radius  =  60  fin  No.  of  particles  =  125 
Particle  radius  =  100  fin  No.  of  particles  =  75 


Run  No. 

'^MD 

(x  10"^^  m^) 

Run  time 

(mins) 

1 

66.034 

212 

2 

56.693 

202 

3 

66.023 

197 

4 

66.043 

210 

5 

65.860 

200 

6 

65.883 

211 

7 

66.063 

205 

Average 

66.9 

205 

Std  Dev 


0.1 
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Table  12  values  for  the  two-particle  size  distribution  simulations 


Radii 

(tm) 

N 

30/50 

125/75 

0.70 

30/100 

125/75 

0.70 

60/100 

125/75 

0.68 

26 


£ri 


Figure  2 


PARTICLE/FLOOR  DISTANCE  (um) 


Force  versus  particle/floor  distance  for  1 00  pm  particles  and  floor  of 
container  interaction,  calculated  using  equation  4. 
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||  30  pm 

\ 

H 


Probability  of  finding  a  particle  as  a  function  of  particle  height  above 
the  floor  of  the  container  for  the  30/50  um  simulation.  The  solid  line 
represents  the  probability  curve  for  the  30  um  particles  and  the  dashed 
line  is  the  probability  curve  for  the  50  ym  particles. 


30 


r 


r 


0.0  175.0  350.0  525.0  700.0 

Height  (jujn) 


Figure  5  Probability  of  finding  a  particle  as  a  function  of  particle  height  above 
the  floor  of  the  container  for  the  60/100  fjtn  simulation.  The  solid  line 
is  for  the  50  iim  particles  and  the  dashed  line  depicts  the 
100  fjm  particles. 
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Probdb i 


Height  (yum) 


Pnbability  of  finding  a  particle  as  a  function  of  particle  height  above 
the  floor  of  the  container  for  the  30/100  um  simulation.  The  solid  line 
and  the  dashed  line  represent  the  probability  curves  for  the  30  and 
100  jjnj  particles  respectively. 
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dynamics  simulation  was  undertaken.  Calculations  were  performed  on  particles  of  30,  50 
and  100  m  radius.  For  particle  size  distributions  where  all  of  the  particles  have  the  same 
radius  (monomodal)  it  was  found  that  the  density  of  the  sediment  is  always  less  than  the 
density  of  an  equivalent  close  packed  structure.  In  addition,  an  increase  in  particle  radius 
leads  to  a  decrease  in  the  sediment  density.  For  bimodal  particle  size  mixes,  either  30/50, 
50/100  or  30/100  m  radii,  fractionation  was  found  to  occur.  In  all  cases,  the  larger 
particles  se'  tied  before  the  smaller  particles  giving  rise  to  a  layer  of  small  particles  on  the 
top  of  the  sediment,  and  a  higher  number  ratio  of  larger  particles  in  the  lower  portions. 
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